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Abstract 

Morphogenesis is a developmental process in which cells organize into 
shapes and patterns. Complex, multi-factorial models are commonly used 
to study morphogenesis. It is difficult to understand the relation between 
the uncertainty in the input and the output of such ‘black-box’ models, 
giving rise to the need for sensitivity analysis tools. In this paper, we 
introduce a workflow for a global sensitivity analysis approach to study 
the impact of single parameters and the interactions between them on 
the output of morphogenesis models. To demonstrate the workflow, we 
used a published, well-studied model of vascular morphogenesis. The 
parameters of the model represent cell properties and behaviors that drive 
the mechanisms of angiogenic sprouting. The global sensitivity analysis 
correctly identified the dominant parameters in the model, consistent with 
previous studies. Additionally, the analysis provides information on the 
relative impact of single parameters and of interactions between them. 
The uncertainty in the output of the model was largely caused by single 
parameters. The parameter interactions, although of low impact, provided 
new insights in the mechanisms of in silica sprouting. Finally, the analysis 
indicated that the model could be reduced by one parameter. We propose 
global sensitivity analysis as an alternative approach to study and validate 
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the mechanisms of morphogenesis. Comparison of the ranking of the 
impact of the model parameters to knowledge derived from experimental 
data and validation of manipulation experiments can help to falsify models 
and to find the operand mechanisms in morphogenesis. The workflow is 
applicable to all ‘black-box’ models, including high-throughput in vitro 
models in which an output measure is affected by a set of experimental 
perturbations. 


Background 


Morphogenesis, the organization of multiple cells into shapes and patterns, is 
a key process in biological development. Computational modeling is commonly 
used to study mechanistic hypotheses on morphogenesis [mnaisisiiiiii as they 
allow for simplification and isolation of the process. These computational studies 
typically involve multi-scale, non-linear and multi-factorial models. So far, the 
behavior of these computational models is studied for one or occasionally two 
parameters at a time, which can lead to a wrong interpretation for non-linear 
models. Studying all parameters collectively with global sensitivity analysis 
resolves this problem. 

In this paper, we introduce a workflow that uses global sensitivity analysis 
to find the relevant single parameters and parameter interactions in ‘black¬ 
box’ models of morphogenesis, which are strongly non-linear and multifactorial. 
Sensitivity analyses of computational models enable us to identify the effects of 
uncertainties in parameter values on the model output. Local sensitivity analy¬ 
sis investigates the behavior of the model in a small region around the nominal 
parameter values and is most often used to study model robustness. Global 
sensitivity analysis (GSA) covers the entire input parameter space, or a specif¬ 
ically selected region hereof. In its most powerful form, it gives information on 
the impact of individual parameters and combinations thereof on a nonlinear 
model for arbitrary parameter distributions. This is what e.g. variance-based 
methods like FAST [5] and the Sobol’ method [HI do. This, of course, is 
computationally expensive, therefore many methods have been proposed with 
simplifying assumptions like linearity of the model (MLR [H]); methods that 
produce less sophisticated results, e.g. partial or no information on interac¬ 
tions (Morris method [12];[Ml[Hj); are less robust like DGSM ([IllUS]); or that 
use prior knowledge of the model, like Bayesian DGSM m)- In this paper 
we use the Sobol’ method cni, where we have modihed the original method 


for efficiency reasons (for more details see Section Global sensitivity analysis). 
Moreover, we introduce an approach to determine the sampling size a priori 
with an a posteriori error check. Thus, it is not likely that the proposed GSA 
will excel in computational efficiency, but it will excel in predictability of the 
costs and reliability of the results. In the biological field GSA is mostly applied 
to models consisting of ordinary differential equations, e.g., in pharmacology 
muni, neurodynamics [20]) or gene expression m and biochemical pathways 
[T7j in cells. 
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GSA can give interesting new insights into models of morphogenesis. Firstly, 
GSA predicts which parameters can best be tuned to affect the model output. 
When the model parameters can be associated with biological cell properties, 
extracellular matrix properties, or gene expression, knowledge of their influence 
on morphogenesis can give predictions for in vitro perturbation experiments, 
e.g. genetic knock-outs. Secondly, apart from identifying the impact of single 
parameters, GSA notably identifies parameter interactions. These can give new 
mechanistic insights in the functioning of the model. Thirdly, GSA is a tool 
to reduce the number of parameters in the model. When the analysis indicates 
that parameter variation does not impact the model output, the parameter value 
can be fixed. Fourthly, GSA can be used to make a selection of models that 
support biologically plausible hypotheses in a set of contradicting mechanistic 
hypotheses. 

As a case study, we performed GSA on a previously published [22], well- 
studied computational model of vascular morphogenesis. In the model, a spheroid 
of cells develops into a vascular network. Gells secrete a compound to which cells 
chemotact by migrating towards higher concentrations of the compound. Vascu¬ 
lar networks form when chemotaxis is inhibited at cell-cell interfaces. Vascular 
endothelial growth factor (VEGF) is a candidate for the secreted compound and 
extensions of pseudopods in the direction of cell-cell contacts might be locally 
inhibited by interference of vascular endothelial cadherins with VEGF receptor 
2 signaling. Because of the key role of such ‘contact-inhibited chemotaxis’ in 
this model, we will henceforth refer to it as the ‘contact inhibition model’. Nu¬ 
merous alternative hypotheses for vascular morphogenesis have been proposed 
ISslEllEaEelEai^EHlElIZS], and it is unsure which of these - if any - is 
correct. Thus the contact inhibition model is here used as an example model 
for morphogenesis, while the proposed GSA approach can assist in falsifying 
mechanisms in the future. 

Figurej^shows the workflow of the GSA analysis proposed in this paper. The 
input (Figure [^), a list of parameter sets, is fed into the cellular Potts-based 
contact inhibition model (Figure]^). This model generates images (Figure]^) 
of the resulting cell configuration as raw output, ranging from spheroids, to 
networks, to dispersed cells. Subsequently, two quantitative output measures 
(compactness and lacuna count) are derived from these images (Figure dP)- 
Two types of GSA are performed on the output measures (Figure [^). Firstly, 
intensity plots show the impact of parameter combinations on the variation in 
the output measures (Figure[^). This analysis only allows for a two-dimensional 
GSA, in which the value of two parameters are varied simultaneously while 
keeping all other parameter values fixed. Secondly, a truly multivariate GSA 
ranks the impact of individual parameters and of parameter combinations on 
the variance of the output measures (Figure[^). Important aspects we address 
in this paper are the reliability and the pitfalls of GSA. 
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G) Mnlti-D analysis of variance 
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individual parameters parameter combinations 

S(p2) = 0.32 S(pi.P4) = 0.1] 

S(p4) = 0.23 S(p2,P3) = 0.06 

S(pi) = 0.11 S(pi,p3) = 0.01 

S(p 3 ) = 0.02 etc. 


Figure 1: Overview of the global sensitivity analysis. (A) The input of 
the model is a list of parameter sets. Each parameter set contains uniformly 
randomly selected values of parameters pi to p^. This input is then fed into (B) 
the Cellular Potts Model (CPM)-based contact inhibition model. (C) The raw 
output of these models are images of the cell configuration at the end of the sim¬ 
ulations. (D) Two output measures, compactness and lacuna count, are derived 
from these images. Two types of global sensitivity analysis are performed on 
these output measures (E). Firstly, intensity plots are used to study the impact 
of two-parameter combinations on the variation in the output measures (F). 
Secondly, Sobol’ indices are used to rank the impact of individual parameters 
and of parameter combinations on the variance of the output measures (G). 
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Methods 


Vascular morphogenesis model 

The contact inhibition model [52] is based on the Cellular Potts Model (CPM) 
|30||3T|. Cells are projected on a regular square lattice (A C as patches of 
connected lattice sites, x. Each lattice site, x G A, that is part of a cell is marked 
with that cell’s identifier (<t(x)) and cell-free lattice sites represent extracellular 
matrix (ECM) with cr = 0. Each cell identifier is associated with a type: T{a) G 
{ECM, cell, border}. Cells have cell properties and behaviors, such as adhesion, 
cell size, or chemotaxis. The forces resulting from the biophysical properties of 
the cells and their active behavior are represented in the Hamiltonian (II) of 
the system 

i7= ^ J(T((7(x)),T(cr(x/))) (1 - S(a(x),a(x/)}} + Aa(o-)^(A(( 7) - a((7))^ , 

(x,x/) cr 

( 1 ) 

in which adhesion (J) is restricted to the cell membrane by the Kronecker delta 
(6(x,y) = {l,x = y;0,x ^ y}) and (a;, x/) represents the set of adjacent lattice 
site pairs. There are three non-zero types of adhesion: Jceii.ceii, -/ceil,border and 
•/cell,ECM- -/cell,cell represents the adhesion strength between cells, and Jceii,ECM 
the adhesion strength of cells to the ECM. The lattice is surrounded with a 
border by which cells are repulsed, by setting Jceii,border = 100. The second 
term constrains the volume of cells, with A representing the preferred size of a 
cell and A/i the rigidity of the cell. Deviation of the actual size (a) of cells from 
their preferred value increases the Hamiltonian. 

A cell moves by copying the state (cr) of a randomly selected lattice site 
X into a randomly selected adjacent lattice site x/. In this manuscript, we 
use the eight, second order neighbors. These copies represent extensions and 
retractions of pseudopods at the cell membrane. A copy attempt that diminishes 
the Hamiltonian represents a move along a force and is always accepted. If a 
copy increases the Hamiltonian, it will only occur due to active movements of the 
cell membrane. We assume that such active motions are distributed according 

-AH 

to a Boltzmann probability function: PBoitzmann(//) = e , with y a measure 
of the amplitude of random membrane fluctuations; /i is a scaling factor of the 
weights (A) of the constraints and we fixed it at /x = 50 conform the settings in 
our previous work [22] • The parameter values can only be qualitatively coupled 
to biological data. 

We assume that cells secrete a chemoattractant at rate a (.s“^), producing 
a concentration field c(x). The chemoattractant diffuses with a diffusion co¬ 
efficient D (mn?/s) and decays with rate e in the ECM, resulting in the 

following partial differential equation (PDE): 

dc 

— = a(l - 5{a{x), 0)) - e5{a{x), 0)c + Z^V^c, (2) 

ot 

such that secretion is located at the cells, where (5(ct(x),0) = 0, and decay in 
the ECM. The field of the chemoattractant is initialized as c(x) = 0 and fixed 
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Figure 2: One-dimensional parameter sweeps for compactness and la¬ 
cuna count. Plots of one-dimensional parameter sweeps for each of the four 
selected parameters: cell rigidity (Ay^), cell-cell adhesion (Jceii.ceii)) diffusion 
coefficient of the chemoattractant (D), and sensitivity of cells to the chemoat¬ 
tractant at cell-matrix interfaces (Ac). The red lines indicate compactness and 
the blue lines lacuna count (mean and standard deviation of 20 simulations). 


boundary conditions are imposed. Cells can respond to this chemoattractant 
by migrating towards higher concentrations (chemotaxis). To this end, the 
change in the Hamiltonian by that copy, AH, is augmented with AiJchemotaxis = 
Ac(c(x) — c(x/)) [35], and contact-inhibition is implemented by setting Ac = 0 at 
cell-cell interfaces such that chemotaxis only occurs at the cell-ECM interface. 

During one time step, referred to as a Monte Carlo Step (MCS), as many 
copy attempts are performed as there are sites in the lattice. The PDF for 
chemoattractant diffusion and degradation (Equation is discretized on the 
CPM lattice and we solve it numerically using a finite-difference scheme. We 
use 15 diffusion steps per MCS. The model is initialized with a centralized 
spheroid of 256 cells within a lattice of 400*400 sites (lattice spacing 2^m). We 
run the model for 5000 MCS, each representing 30 seconds, as networks are well 
formed in this time in the model as well as in vitro [25] . 
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Figure 3: Two-dimensional intensity plots of lacuna count. The intensity 
of the output measure lacuna count, mapped to an interval of 0 to 15 as indi¬ 
cated by the color bars, is plotted for each two-parameter combination of the 
parameters cell rigidity (A^), cell-cell adhesion (Jceii.ceii), diffusion coefficient 
of the chemoattractant (D), and sensitivity of cells to the chemoattractant at 
cell-matrix interfaces (Ac). 
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Figure 4: Two-dimensional intensity plots of compactness. The intensity 
of the output measure compactness, mapped to an interval of 0 to 1 as indi¬ 
cated by the color bars, is plotted for each two-parameter combination of the 
parameters cell rigidity (Aa), cell-cell adhesion (Jceii.ceii), diffusion coefficient 
of the chemoattractant {D), and sensitivity of cells to the chemoattractant at 
cell-matrix interfaces (Ac). 




















Global sensitivity analysis 

The variation in a solution or a measure thereof, like compactness, over the 
complete parameter space can only be visually inspected by looking at one or 
at most two parameters at a time while keeping the others fixed (cf. Figure]^ 
and 1^. Since different measures will produce different multivariate output 
distributions and therefore also might result in a different outcome of the GSA, 
it is important to choose a measure or, more likely a number of measures that 
are significant for the study at hand. If one wants to take the influence of all 
parameters simultaneously into account some form of a global measure of the 
multivariate output distribution is required. One such a measure is the variance 
of the distribution, which will be used in this paper. 

We are specifically interested whether parameter interactions have a large 
impact on the output of this specihc CPM-based model. Interactions of the 
parameters are unpredictable in non-linear models such as the CPM, but their 
impact is significant, since a large combined effect of parameters on the output 
impedes the experimental testing of a predicted effect of a single parameter. 

Sobol’ m Uni introduced so-called global sensitivity indices that describe the 
impact of specific parameters or combinations thereof on the uncertainty in the 
model output and more in particular on the variance of the output distribution, 
hence the term ‘variance-based’ GSA. In the original method the necessary in¬ 
tegrals are computed with Monte Carlo. However, the Sobol’ indices can be 
computed very efficiently when the distribution of the output measure or re¬ 
sponse surface is expanded into a series of orthogonal polynomials, the Polyno¬ 
mial Chaos Expansion (PCE) [33l [34] . An overview of this method, using the 
same notation as in the following, can be found in |35j ; here we summarize only 
the most important definitions for the case that the stochastic input consists of 
independently uniformly distributed random variables. Note, however, that the 
method can also be applied for arbitrary, even non-parametric, distributions, 
allowing for data-driven GSA (see also [36]). The strength of the method de¬ 
scribed below is that it (i) can efficiently study multiple output measures derived 
from the output images, (ii) can robustly identify parameter interactions, and 
(iii) checks the reliability of the result. 

Let ^ be the n-dimensional vector of the independently uniformly distributed 
input parameters and g{^) its joint probability density function (pdf). The 
output measure u{^), e.g., of the (black-box) Cellular Potts Model, is expanded 
into a truncated series of polynomials that are orthogonal with respect to the 
pdf p, separating the output into a deterministic and a stochastic part 

N 

i=0 


where the n-variate polynomials $^(4) are products of n univariate Legendre 
polynomials. The number of expansion terms N is given by -I- 1 = ; 

with n the number of parameters and the approximation order p the highest 
order of <I>i. 
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To compute the expansion coefficients Ui of Equation ([^ we apply Spectral 
Projection which has the advantage that it can be used for black-box models 
since it projects the solution - and not the model - onto the polynomial space 

where 5 is the support of the joint pdf g{^). As the parameter inputs are 
independent, both and g can be written in product form; for the multivariate 
polynomial this results in a product of univariate polynomials 

n 

®i(€) = n ®index(*,fc)(Cfe), with index(i, fc) = {0, ...,p} and ^'o(Cfe) = 1- (5) 

k=l 

The integrals in Equation Q can then be computed by a repeated one-dimensional 
Gauss-Legendre quadrature rule 

^ W, Af, n 

~ 11 112 'y ' ' ’ ’ ^ ' ^(Czi J ■ ■ ' I ’^Ik '^*index(z,fe) (Cu )^ (6) 

*11 ii = l l„ = l fc=l 


with Nq the number of quadrature points and w the associated weights. Note, 
that for integrals with a known weight function, like e.g. a pdf. Gauss quadrature 
has the optimal convergence order of 2 Nq — 1 for Nq quadrature points, where 
the points and the weights of the quadrature rule are dependent on the weight 
function. 

How to choose Nq and p to obtain reliable Sobol’ indices will be the subject 


of Section Reliable GSA in practice 


Statistics and Polynomial Chaos Expansion 


Using a PC expansion, the only input needed to compute the moments and the 
Sobol’ indices of the output distribution are the expansion coefficients. E.g., 
the mean g, = uq and the variance is given by 


N 

J (u(0 - gf gi^)d^ « =: VarpcE 


(7) 


Note, that the approximation, VarpcE, is a monotonously increasing function 
of N and thus of p. The sum in the variance formula can be directly split into 
contributions from the various parameters or combinations thereof, the Sobol’ 
indices (cf. [37]). E.g., for the first-order Sobol’ index for parameter j only 
terms contribute if $^(4) equals a univariate polynomial in 


E^lbool(bj) 

VarpcE 


5, 


( 8 ) 


where bool(z,j) = (index(z, j) > 0 A index(i, fc) = 0,Vfc ^ j). For a combined 
influence of more than one parameter like S 13 the Sobol’ index can be computed 
analogously. The sum of all Sobol’ indices equals one. 
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Reliable GSA in practice 

At first sight the accuracy of the PCE approximation of the response surface - 
and thus of the statistics - seems to be determined by the number of expansion 
terms, N, in Equation (|^. But the accuracy of the expansion coefficients ut 
also plays an important role. This accuracy is determined by the approximation 
(Equation]^ of the integral in Equation Q, which is determined by the number 
of quadrature points, Nq. Moreover, the higher PCE order p needed to obtain 
sufficient accuracy, the higher the polynomial order of becomes, which 

increases the complexity of the integrand. If one computes the integral of a 
high order polynomial with a small amount of points, the resulting expansion 
coefficients are merely noise instead of being informative. 

The question we want to answer in this section is how to determine the 
number of quadrature points, Nq, and the expansion order, p, to obtain a suffi¬ 
ciently high accuracy for the coefficients to allow us to trust the Sobol’ indices 
and more specifically the ranking of the parameters that follows from it. Here, 
we sketch a method to determine the number of quadrature points. It relies on 
the fact that the Sobol’ indices are variance-based, i.e., one can not expect to 
compute Sobol’ indices accurately from a PCE expansion for which the variance 
(Equation is not a sufficiently accurate approximation of the true value or 
at least comparable to the Gauss-Legendre quadrature approximation of the 
integral. So, let us define 


errvar := Vardata - VarpcE, 


(9) 


Vardata — 

E" 

• E • 


Mdata — 

h=i 

N, 

E" 

• E • 



h=i 




For a given choice of Nq one can easily compute PC expansions for various 
orders p. If erryar is small and the required Sobol’ indices have converged, the 
result can be trusted. 

We illustrate this approach with a function for which the values of the statis¬ 
tics are analytically known, viz., the Ishigami function [381139) 

/(I) = sin(^i) -h asin^(^ 2 ) + &C 3 sin(Ci), (10) 

with ~ TT, tt], i = {1, 2, 3}, and a = 7 and b = 0.1. We compute the PCE 
approximation of Equation ( |10[ ) for an increasing number of PCE terms and an 
increasing number of quadrature points. Table illustrates the result of using 
not enough quadrature points {Nq = 2 and Nq = 5): there is no convergence 
in the statistics of the PCE approximation and for p = 3 and 6, respectively, 
the noise takes over and the results are meaningless. Tableshows that, using 
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Table 1: Statistics computed with insufScient quadrature points. The resulting 


PCE approximation and thus the statistics can not be trusted. 



P 

V^^data 

VarpcE 

^1 

^2 

Sl 3 

^3 

Si2 

S23 

2 

1 

4.09 

4.09 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 


2 

4.09 

4.09 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 


3 

4.09 

8.32 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 


4 

4.09 

185.91 

0.36 

0.32 

0.00 

0.32 

0.00 

0.00 


5 

4.09 

197.44 

0.34 

0.30 

0.03 

0.30 

0.03 

0.00 

5 

1 

18.60 

2.64 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 


2 

18.60 

2.70 

0.98 

0.02 

0.00 

0.00 

0.00 

0.00 


3 

18.60 

6.22 

0.69 

0.01 

0.30 

0.00 

0.00 

0.00 


4 

18.60 

17.17 

0.25 

0.64 

0.11 

0.00 

0.00 

0.00 


5 

18.60 

18.50 

0.23 

0.60 

0.17 

0.00 

0.00 

0.00 


6 

18.60 

29.49 

0.14 

0.75 

0.11 

0.00 

0.00 

0.00 


7 

18.60 

31.41 

0.19 

0.70 

0.11 

0.00 

0.00 

0.00 


8 

18.60 

31.50 

0.19 

0.70 

0.11 

0.00 

0.00 

0.00 


9 

18.60 

37.51 

0.23 

0.59 

0.18 

0.00 

0.00 

0.00 


10 

18.60 

88.82 

0.29 

0.43 

0.08 

0.20 

0.00 

0.00 

Exact 

13.84 

0.31 

0.44 

0.24 

0 

0 

0 


sufficient quadrature points, for an increasing number of expansion terms the 
PCE variance converges to the data variance. If both variances are alike also the 
Sobol’ indices have converged to the true values. Still the number of expansion 
terms should not be taken too large as can be seen for p > 9 and p > 13 where 
again the noise gradually takes over. 

Finally, we also used the Saltelli method[40] - an improvement of the original 
Sobol’ method - to compute the Sobol’ indices for this problem. To reach a sim¬ 
ilar accuracy approximately 70-80 times as many sampling points are required, 
thus showing the gain in efficiency using the PCE-Gauss method to compute 
the Sobol’ indices. 

Software and computational dataset 

All software used in this paper is publicly available. The contact inhibition 
model resides at http://sourceforge.net/projects/tst/. For GSA we pro¬ 
vide a repository containing the computational dataset and the analysis soft¬ 
ware at http://persistent-identifier.org/?identifier=urn:nbn:nl:ui: 
18-23590, 

Results 

As a case study for the global sensitivity analysis (GSA) approach, we used a 
well-studied computational model of vascular morphogenesis: the contact inhibi¬ 
tion model |22j . We studied what single parameters and parameter interactions 
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Table 2: Statistics computed with sufficient quadrature points. The PCE ap¬ 


proximation and the statistics show convergence (bold lines). 


Nq P 


VarpcE 

Si 

^2 

5i3 

^3 

Si2 

S 23 

8 1 

13.59 

2.64 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 

2 

13.59 

3.00 

0.88 

0.12 

0.00 

0.00 

0.00 

0.00 

3 

13.59 

6.55 

0.66 

0.05 

0.29 

0.00 

0.00 

0.00 

4 

13.59 

10.34 

0.42 

0.40 

0.18 

0.00 

0.00 

0.00 

5 

13.59 

11.73 

0.37 

0.35 

0.28 

0.00 

0.00 

0.00 

6 

13.59 

13.45 

0.32 

0.44 

0.24 

0.00 

0.00 

0.00 

7 

13.59 

13.59 

0.32 

0.43 

0.25 

0.00 

0.00 

0.00 

8 

13.59 

13.59 

0.32 

0.43 

0.25 

0.00 

0.00 

0.00 

9 

13.59 

13.59 

0.32 

0.43 

0.25 

0.00 

0.00 

0.00 

10 

13.59 

15.32 

0.28 

0.50 

0.22 

0.00 

0.00 

0.00 

11 

13.59 

15.35 

0.29 

0.49 

0.22 

0.00 

0.00 

0.00 

12 

13.59 

19.17 

0.23 

0.60 

0.18 

0.00 

0.00 

0.00 

13 

13.59 

21.08 

0.29 

0.54 

0.17 

0.00 

0.00 

0.00 

14 

13.59 

21.51 

0.28 

0.55 

0.17 

0.00 

0.00 

0.00 

15 

13.59 

27.62 

0.32 

0.43 

0.25 

0.00 

0.00 

0.00 

10 1 

13.84 

2.64 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 

2 

13.84 

3.00 

0.88 

0.12 

0.00 

0.00 

0.00 

0.00 

3 

13.84 

6.54 

0.66 

0.05 

0.29 

0.00 

0.00 

0.00 

4 

13.84 

10.36 

0.42 

0.40 

0.18 

0.00 

0.00 

0.00 

5 

13.84 

11.75 

0.37 

0.35 

0.28 

0.00 

0.00 

0.00 

6 

13.84 

13.59 

0.32 

0.44 

0.24 

0.00 

0.00 

0.00 

7 

13.84 

13.72 

0.32 

0.44 

0.25 

0.00 

0.00 

0.00 

8 

13.84 

13.84 

0.31 

0.44 

0.24 

0.00 

0.00 

0.00 

9 

13.84 

13.84 

0.31 

0.44 

0.24 

0.00 

0.00 

0.00 

10 

13.84 

13.84 

0.31 

0.44 

0.24 

0.00 

0.00 

0.00 

11 

13.84 

13.84 

0.31 

0.44 

0.24 

0.00 

0.00 

0.00 

12 

13.84 

13.95 

0.31 

0.45 

0.24 

0.00 

0.00 

0.00 

13 

13.84 

13.95 

0.31 

0.45 

0.24 

0.00 

0.00 

0.00 

14 

13.84 

15.81 

0.27 

0.51 

0.21 

0.00 

0.00 

0.00 

15 

13.84 

15.84 

0.28 

0.51 

0.21 

0.00 

0.00 

0.00 

Exact 

13.84 

0.31 

0.44 

0.24 

0 

0 

0 
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are important in the development of a spheroid of cells into vascular networks. 
For this purpose, we used the procedure outlined in Figure 1) select output 
measures, 2) select input parameters, 3) select a relevant subset of the global 
parameter space, 4) analyze the raw output, 5) perform GSA. 


Selection of output measures 

The contact inhibition model [55] produces images of cell configurations as raw 
output. We chose two measures to quantify the raw output: compactness and 
lacuna count. Compactness of the network is a suitable measure of network 
development [55] and is defined as the ratio Aceiis/Ahull, with Aceiis the number 
of lattice sites occupied by cells within a convex hull around all cells and Ahull 
the total number of lattice sites within the convex hull. A solid spheroid and 
a confluent monolayer of cells have a compactness close to one, while networks 
that contain lacunae have low values for compactness. Lacuna count is the 
number of lacunae in a network. Lacunae are defined as patches of medium 
(connected components of cr = 0) enclosed by cells and are only counted when 
they have at least the size of a cell (50 lattice sites ~ 200 


Selection of input parameters 


The contact inhibition model [55] is a stochastic, multi-factorial model. We 
refer to Section Methods for a detailed description of the model. The contact 


inhibition model has nine parameters: the number of cells (N)^ the target size 
of a cell (A), the rigidity of the cell (A^), cell-cell adhesion (Jceii.ceii), adhe¬ 
sion between cells and the extracellular matrix (Jceii.ECM), the secretion rate 
of a chemoattractant by cells (a), a diffusion coefficient of the chemoattractant 
(D), the decay rate of the chemoattractant (e), and a sensitivity of cells to the 
chemoattractant at cell-matrix interfaces (Ac). 

In total, there are four model components or mechanisms in the contact in¬ 
hibition model, namely cell size, adhesion, contact-inhibited chemotaxis and the 
gradient of the chemoattractant. In order to study the impact of each mech¬ 
anism in the model extensively, we selected one parameter for each, ensuring 
that it is computationally feasible to generate enough data points for reliable 
GSA results. We thus selected four parameters: the cell rigidity (Aa), cell-cell 
adhesion ( Jceii.ceii), the diffusion coefficient of the chemoattractant (D), and a 
sensitivity of cells to the chemoattractant at cell-matrix interfaces (Ac). The 
other parameters that regulate cell size (A), adhesion (Jceii.ECM), or the gra¬ 
dient of the chemoattractant (e and a) will be fixed at the reference values 
corresponding to the values in [55] . We kept the number of cells {N) in the 
spheroid constant, because we know from experience that it does not influence 
sprouting of spheroids in our model. 

A GSA with four parameters can give new insights as four parameters are too 
many to obtain the relative impact of the parameters and their interactions with 
visual plots or to know their effect solely by logic or intuition, while the number 
of simulations required for a GSA with four input parameters is computationally 
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Table 3: Overview of the parameter selection for the GSA. The names of the 
parameters are listed in the first column, a parameter description in the second 
column, and the selected parameter value ranges in the last column. 


Name 

Description 

Range 

Ac 

Chemotaxis 

10 to 3000 

D 

Diffusion coefficient 

le-14 to 5e-13 


A Area 

5 to 300 

^ce\l,ce\l 

Cell-cell adhesion 

0 to 120 


very feasible. A GSA with all parameters of the model is not expected to give 
additional information on the relative balance of the mechanisms and would 
be very time-consuming for a computationally intense model like the contact 
inhibition model. It would require roughly 10® simulations (c.f.. Section Reliable 


GSA in practice) to obtain reliable GSA results with all model parameters. 


Selection of a relevant subset of the global parameter space 

To select the parameter ranges for which spheroids of cells develop into net¬ 
works, we studied one-dimensional parameter sweeps of the four selected input 
parameters for the compactness and lacuna count (Figure]^. The red lines in 
Figure represent the compactness and the blue lines the lacuna count. We 
selected the region in which the morphology of the network, and thus the value 
of the output measures, is changing and where no model artefacts arise. It is 
well studied for which parameter ranges artefacts arise in the CPM [B], such 
as lattice anisotropy and ‘frozen’ motility of cells. The regions shaded in gray 
indicate the deleted regions from the parameter space. For Xa the region 0 
to 5 is deleted: cells cannot retain their volume here and disappear. This is 
a model artefact and does not represent a biological plausible situation. The 
region Aa > 300 is deleted, because cells are so rigid here that they hardly 
move. In the region Ac < 10 only spheroids form and for Ac > 3000 similar 
networks are always formed, thus these regions are deleted because the network 
morphology does not change. The parameters and their selected value ranges 
are listed in Table [H 


Based on the reliability study for the Ishigami test model (see Section Re¬ 
liable GSA in practice), we expected that we required 10000 data points to 
perform a reliable GSA on our model. The points were chosen according to 


the Gauss-Legendre quadrature rule (see Section Global sensitivity analysis), 
resulting in ten values for each parameter. To correct for the stochasticity of 
the contact inhibition model, each parameter set is replicated twenty times with 
a different random seed and the output is averaged over them. The size of the 
standard deviations in Figure indicate that the variation over different ran¬ 
dom seeds is very small for compactness, whereas the stochasticity in the model 
has a larger affect on the lacuna count. Nevertheless, this lacuna count is a 
reasonable measure for the network morphology. 
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Figure 5: Overview of the raw output. An overview of the raw output 
of the contact inhibition model, the cell configurations at the end of a simu¬ 
lation, is shown in a collage of images. The cell rigidity (Aa) is varied over 
the horizontal axis and cell-cell adhesion (Jceii.ceii) over the vertical axis. For 
each selected combination hereof, a subcollage is shown in which the diffusion 
coefficient of the chemoattractant {D) is plotted against the sensitivity of cells 
to the chemoattractant at cell-matrix interfaces (Ac). 
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Analysis of the raw output 

The raw output of a model simulation is an image of the cell configuration 
at the end of a simulation. Figure gives an overview of the raw output for 
the selected parameter space. Examples of possible morphologies are shown in 
Figure ranging from spheroids to small networks with one lacuna, to fine- 
mazed networks with many lacunae. This is a visual reassurance that the input 
parameter space is well chosen. However, it is very difficult to predict from 
the raw output which parameters have a strong impact on the development of 
networks from spheroids. A GSA can give us insights into this, as we will show 
in the next section. 


GSA of uetwork developuieut froui spheroids 

We performed two types of GSA on the distribution of the output measures, 
compactness and lacuna count, to study the impact of the parameters on vas¬ 
cular network development. The first type of GSA studies the variation of the 
output measures and the second type studies the decomposition of the variance 
of the distribution of the output measures. 

GSA of the variation of the output measures 

The variation in the output measures can be visualized by plotting the inten¬ 
sity of the output measures over two-dimensional slices of the parameter space. 
Figure and show the intensity plots of the lacuna count and compactness, 
respectively, for each possible pairing of parameters. The parameter values are 
selected according to the Gauss-Legendre quadrature rule. 

Figure shows that the diffusion coefficient is the main source of variation 
for the lacuna count: the lacuna count is high for low values of D and low 
for high values of D, independent of the other parameters. That the lacuna 
count does not vary significantly over the entire perpendicular axis indicates 
that the parameter of the perpendicular axis does not have much impact. The 
dominance of the diffusion coefficient masks the impact of the other parameters. 
To reveal the impact of the other parameters. Additional file 4: Figure SI caps 
the intensity values at a lacuna count of five. 

Figure shows a high variation of the compactness in each plot. As a 
consequence, it is difficult to determine which parameters have a dominant 
impact on compactness. Interactions between parameters are difficult to deduce 
from these two-dimensional intensity plots. A variance-based GSA is well suited 
to derive parameter interactions and the ranking of individual parameter effects, 
as will be outlined in the following subsection. 

Variance-based GSA of the output measures 

To study the impact of single parameters and of parameter combinations on the 
development of networks from spheroids, we performed a GSA of the output 
distribution of compactness and lacunae count using the Sobol’ indices. We 


17 


refer to Section Global sensitivity analyst^ for a detailed description of how to 
obtain the Sobol’ indices that represent the impact of the parameters. The GSA 
results of both measures are reliable, since the Sobol’ indices have converged for 
values of p for which erryar (Equation]^ is small (see Additional file 1: Table 
SI and Additional file 2: Table S2). 


Table 4: Global sensitivity analysis results. The Sobol’ indices for the individual 
parameters (indices above mid-line) and for their combinations (indices below 
mid-line) are listed for the GSA of compactness and lacuna count. 



compactness 

lacuna count 

S(Ac) 

0.3188 

0.0074 

S{D) 

0.2969 

0.7130 

S{Xa) 

0.0266 

0.0125 

S(</cell,cell) 

0.2048 

0.0407 

S(Ac,Z?) 

0.0125 

0.0570 

S(Ac,Aa) 

0.0107 

0.0043 

S ( Ac, J cell,cell) 

0.0559 

0.0145 

SiDM) 

0.0017 

0.0347 

S (Z?, (/cell,cell) 

0.0127 

0.0521 

S( A^,t7cell,cell) 

0.0102 

0.0048 

S(Ac,Z?,Aa) 

0.0102 

0.0315 

S(Ac,/^,'/cell,cell) 

0.0257 

0.0232 

S ( Ac, Ay4, (/cell, cell) 

0.0217 

0.0131 

S (/^, Ayl, J cell,cell) 

0.0075 

0.0094 


The second column of Table lists the impact of the individual parameters 
and their combinations on compactness. The sensitivity for the chemoattractant 
at cell-matrix interfaces (Ac) has the highest impact on network development 
(S(Ac)=0.3188), followed by the diffusion coefficient with S(iA)=0.2969, and 
cell-cell adhesion with S( Jceii,ceii)=0.2048. Elasticity of cells has a low impact of 
(S(Aa)=0.0266). Seventeen percent of the variance is caused by interactions of 
parameters. Ac and Jceii.ceii have a combined impact of 0.0559. The impact of all 
other interactions was lower than S(Aa), which we will consider as a threshold 
for relevant impact. 

The third column of Table lists the impact of the individual parameters 
and their combinations on lacuna count. The individual impact of the diffu¬ 
sion coefficient is dominant, with S(D)=0.7130. Gell adhesion also has a small 
individual impact (S(JceU,ceii)=0.0407). In total, twenty four percent of the 
variance is induced by combinations of parameters. There are five parameter 
combinations, which all include the diffusion coefficient, with a higher impact 
than the threshold: S(Ac,T>)=0.0570, S(Z?,Aa)= 0.0347, S(i:),JceU.ceii)=0.0521, 
S(Ac,7I, Jcen,ceii)=0.0476, and S(Ac,7?,A^)=0.0315. The total impact of the dif¬ 
fusion coefficient is 90 percent. When we focus on low values of the lacuna count, 
by capping the lacuna count at a maximum of five lacunae, the dominance of 
the diffusion coefficient is slightly reduced and an extra interaction of Ac and 
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Figure 6: Dominant effect of the diffusion coefficient on lacuna count. 

A collage of the cell configurations at the end of simulations in the contact 
inhibition model, in which the diffusion coefficient of the chemoattractant (D) is 
varied over the horizontal axis and the sensitivity of cells to the chemoattractant 
at cell-matrix interfaces (Ac) over the vertical axis. 


'^cell.cell is found (Ac,dcell,cell) —0- 0409) (Additional file 3: Table S3). 


Interpretation of the GSA results 

The GSA results show that three parameters account for over 80 percent of the 
variance of the compactness distribution. Consistent with previous studies of the 
contact inhibition model [12], these three parameters are the diffusion coefficient, 
sensitivity to the chemoattractant at cell-matrix interfaces and cell-cell adhesion. 
For the lacuna count, the GSA identified solely the diffusion coefficient as the 
dominant parameter. This dominant effect is apparent in a collage of output 
images (Figure]^: the number of lacunae varies over the horizontal axis that 
represents the diffusion coefficient D, whereas there is little variation along 
the vertical axis that represents the sensitivity to the chemoattractant at cell- 
matrix interfaces Ac. The number of lacunae is the largest for small values of 
the diffusion coefficient (around D = 4.3* 10“^^ rn?/s), whereas no lacunae are 
formed for large values of the diffusion coefficient. A similar trend is seen when 
the diffusion coefficient is plotted against cell-cell adhesion or cell rigidity (not 
shown). The distance over which adjacent branches can attract one another 
is given by the length of the chemoattractant gradient (Equation , which is 
characterized by the diffusion length, Ld = yJTf~D^ the distance over which the 
secreted chemoattractant drops to 1/e of the concentration at the cells (see, e.g., 
the discussion in Ref. [22|)- If Ad becomes shorter, branches that would fuse for 
larger values of Ad will not fuse. Hence the pattern will be more fine-grained. 
Also a shorter value of Ad will create sharper gradients and hence increase 
the inward chemotactic force (as AiJ = Ac * gradient) hence ’’squeezing” the 
branches more and making them thinner. 


19 




In conclusion, the GSA is able to identify the dominant single parameters 
for compactness and lacuna count. In addition, it gives new information on the 
relative ranking of the impact of these single parameters. 

In contrast to the one-dimensional parameter studies performed in |22j , GSA 
provides information on interactions of parameters. Combinations of parameters 
account for 17% of the variance in the compactness distribution, and for 24% 
of the lacuna count distribution. This indicates that most parameters impact 
the model output independently. Interestingly, the parameter combination of 
Ac and Jceii.ceii impacted the lacuna count (as seen in Additional hie 4: Figure 
SI, which is capped at a maximum of hve lacunae) as well as compactness. How 
can we explain this interaction? Sprout formation requires a balance between 
Ac-dependent chemotaxis, creating an inward force perpendicular to the sprout 
surface, and Jceii,cell-dependent cell-cell adhesion, which is responsible for the 
surface tension of individual cells. In the limit of zero-surface tension, the 
cells would behave as a zero-viscosity huid and the chemotaxis would compress 
sprouts until they become inhnitely thin |41j . The cellular surface tensions resist 
such compression, thus determining the thickness of sprouts. Altogether, this 
parameter interdependence highlights a new insight in the mechanisms driving 
sprouting in our model. 


Discussion 

Biological morphogenesis is a highly complicated process, involving genetic 
regulation, pattern formation, the biophysics of collective cell migration, me¬ 
chanical cell-cell interactions, and so forth. As such multiscale mechanisms 
are practically impossible to understand intuitively, in recent years model¬ 
ing and simulation has become a key tool in developmental biology (see, e.g., 
refs. 132132031 US]). These efforts have led to highly complicated models, where 
traditional analysis tools in dynamical systems theory, such as bifurcation anal¬ 
ysis and phase plane analysis, fall short. The models must then be treated as 
‘black-box’ systems: one- or two-dimensional parameter sweeps are performed, 
creating images and movies as output, which can be used to obtain various 
quantitative output measures. These parameter sweeps must be started from 
one or a few nominal parameter sets around which n-dimensional cross-shaped 
sweeps through the parameter space are performed. However insightful such 
studies are, a danger is that the impact of some parameters is overlooked: the 
conclusions may depend on what sets of nominal parameter values were selected. 
Using a simple, published simulation model of vascular morphogenesis, we have 
shown in this work how a multivariate GSA helps to get more insight in the 
relative impact of single parameters and of their interactions. We introduced a 
workflow for GSA of ‘black-box’ models of morphogenesis. 

We applied the workflow to a vascular morphogenesis model which we refer 
to as the ‘contact inhibition model’. The output of the contact inhibition model 
consists of images of the cell configuration in a simulation. To quantify network 
development, we measured the compactness and the lacuna count of the cell 
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configuration at the end of the simulation. A GSA with four input parameter 
distributions, that each described one of the four general model components, 
was performed for both measures. The GSA results of compactness and lacuna 
count both indicated that variation of the rigidity of the cells (A^i) has very 
little impact on the model output. As a result, the model can be reduced by 
fixing this parameter. For compactness, the sensitivity for the chemoattractant 
at cell-matrix interfaces (Ac) has the highest impact on network development, 
followed by the diffusion coefficient (D) and cell-cell adhesion (Jceii.ceii)- In con¬ 
trast, the GSA showed that the diffusion coefficient alone is dominant for lacuna 
count. The results for both measures are in line with what has been previously 
reported [H]. New information from the GSA results is the relative impact 
of the single parameters. In addition, GSA identified interactions between pa¬ 
rameters, which have led to new insights in the mechanism of sprouting in the 
model. Most notably, the parameter interactions in this specific GPM-based 
model have very low impact. Since GSA has not been performed for GPM- 
based models before, it is an important new insight for the GPM community 
that the most basic mechanisms of the GPM, such as cell size and adhesion here 
function independently. 

Besides the contact inhibition model, there are multiple other computational 
models of vascular network development [531 HH113 HSl HZl HE] ■ These models 
often share common mechanisms that drive sprouting, but differ by one or a few 
mechanisms. It is still not known which mechanisms drive sprouting in vivo^ 
or whether a different set of mechanisms is used under different conditions. 
We propose GSA as an approach to help falsify these models. Firstly, the 
ranking of the relevance of the mechanisms in the models can be compared 
with knowledge of the impact of these mechanisms from experimental data to 
falsify models. A second model falsihcation method is the validation of the 
experimental predictions of each model based on the GSA results. 

The workflow is designed to take into account some pitfalls of GSA. These 
arise from the dependency of the outcome on the choices one makes for the 
output measure, input parameters and their distributions. Different output 
measures can give different results, as was the case for compactness and lacuna 
count. This indicates that it is essential to consider carefully whether the se¬ 
lected measure truly describes your goal and if there are other measures for it. 
A selection of input parameters might be necessary when it is not computation¬ 
ally feasible or methodologically desirable to use all parameters of the model. 
The importance of the selection of the correct parameter distributions has also 
been discussed elsewhere m- Intuitively, a large range for the parameter values 
will allow for the largest variation in the output and thus the most interesting 
result. However, since the analysis is global over the entire parameter space, 
local though important features might become unnoticed if the distribution is 
too widespread. For instance, for the contact inhibition model we were inter¬ 
ested in the region where the networks developed and where the measures were 
changing accordingly, and variation in these regions could become unnoticed if 
we included large regions where for instance spheroids do not sprout. Ideally, 
the parameter distribution comes from experimental measurements, but in ab- 
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sence hereof we propose to study the variation of the output measures for each 
parameter individually. 

It is crucial to have an estimate of the accuracy of the sensitivity results. One 
option is to compare the results with the outcome of an analysis with a higher 
accuracy computed with more quadrature points and a higher PCE order, like 
advocated in m- In this paper we proposed a simpler and cheaper rule: given 
the number of quadrature points the Sobol’ indices should show convergence 
for those values of p for which the variance computed with the Gauss-Legendre 
quadrature rule is more or less equal to the variance computed from the PCE 
approximation. If a higher PCE order is required, more model simulations are 
needed. Since the computation of the PCE statistics is ‘for free’ compared to 
model simulations this is an efficient way of judging whether the accuracy of 
the statistics is sufficient for one’s aim. Although Gauss quadrature is optimal, 
it has the disadvantage that it is not a nested quadrature rule, i.e., if more 
quadrature points are required, the old model results cannot be re-used. An 
alternative for Gauss quadrature is Monte Carlo (MC) integration. Sampling 
the PCE integrals by MC is less optimal, so more simulations are needed to 
obtain reliable GSA results. For the Ishigami test model, MC needs a 100 times 
more simulations to get comparable results. The benefit of MC is that you 
can check ‘on the fly’ if there are enough data points generated to get reliable 
results. Adding simulations on the fly is particular useful when the estimated 
number of simulations based on the Gauss quadrature rule is computationally 
unfeasible, but one expects or hopes that the output distribution is relatively 
smooth and thus can be described by a low order PCE approximation. 

Some studies require GSA of a subspace of the output distribution. For 
instance in our case study, to study not the conditions for network formation 
per se, but the details of the network morphology (e.g. branch length, branch 
thickness, and so forth), we must preselect a region of the parameter space 
where networks actually form. Unfortunately, such a subspace would no longer 
guarantee that the input distribution is independently random uniform. For 
such cases, a more complicated method to compute the Sobol’ indices [3S] is 
required. 

Besides in computational models, the impact of biological factors on mor¬ 
phogenesis is also studied in vitro. High-throughput image-based screenings 
systematically analyze the impact of genes or potential drugs on cell behavior, 
such as cell migration |46j . This ‘systems microscopy’ approach is well suited for 
parallel screening of cellular responses to numerous experimental perturbations 
m- Such high-throughput screens can be performed for the genes, growth fac¬ 
tors or ECM concentrations affecting morphogenesis. This is conceptually very 
similar to parameter studies of in silico ‘black-box’ models. The perturbed bio¬ 
logical factors represent the input parameters and the output is an image from 
which quantitative data can be derived. Therefore, the GSA workflow proposed 
in this paper is directly applicable to high-throughput in vitro studies. 
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Conclusions 


Morphogenesis is a complex biological process in which cells organize into shapes 
and patterns. Computational modeling is used to get insights in the mecha¬ 
nisms of morphogenesis. These models are often multi-scale, non-linear and 
multi-factorial, making it difficult to relate their input to their output. The 
behavior of such ‘black-box’ models is mostly studied by visual inspection and 
analyses of the individual output (e.g. images and movies) and with one- or 
two-dimensional parameter sweeps of output measures. However, this does not 
provide insight in the relative impact of single parameters and of their interac¬ 
tions on the outcome of the model. GSA fulhlls this task. GSA results can give 
insights in the dynamics of the model and help to generate experimental pre¬ 
dictions to manipulate morphogenesis. In this paper, we introduced a workflow 
for GSA of such models and addressed pitfalls and reliability of the analysis. 
The workflow is applied to the contact inhibition model, a cell-based model of 
vascular morphogenesis. GSA was able to correctly identify dominant param¬ 
eters and gave new insights on the magnitude and ranking of their individual 
impact and importantly, on their interactions. In summary, we propose GSA of 
‘black-box’ models, such as complex computational models or high-throughput 
in vitro models, as an alternative approach to get insights in the mechanisms of 
morphogenesis. 
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Supplementary Information 

Additional file 1 — Global sensitivity analysis results for 
compactness 

Table SI. Global sensitivity analysis results for compactness. 


p 

12 

13 

14 

15 

Variance data 

0.0453 

0.0453 

0.0453 

0.0453 

Variance PGE 

0.0452 

0.0453 

0.0454 

0.0456 

S(Ae) 

0.3190 

0.3188 

0.3186 

0.3183 

SiD) 

0.2971 

0.2969 

0.2965 

0.2958 

S(Aa) 

0.0267 

0.0266 

0.0266 

0.0265 

S(</cell,cell) 

0.2052 

0.2048 

0.2043 

0.2032 

S{Xc.D) 

0.0124 

0.0125 

0.0126 

0.0127 

S{Xc.Xa) 

0.0107 

0.0107 

0.0109 

0.0110 

S(Ac5t/cell,cell) 

0.0558 

0.0559 

0.0561 

0.0570 

S{DM) 

0.0016 

0.0017 

0.0017 

0.0017 

S(Z?, (/cell,cell) 

0.0127 

0.0127 

0.0127 

0.0126 

S( A^,t/cell,cell) 

0.0102 

0.0102 

0.0102 

0.0101 

S(Ac,i/,A^) 

0.0098 

0.0102 

0.0104 

0.0107 

S ( Ac, /^, (/cell,cell) 

0.0253 

0.0257 

0.0262 

0.0268 

S ( Ac, Aa^ ,t/cell,cell) 

0.0214 

0.0217 

0.0220 

0.0225 

S (-Z/, A , (/cell, cell) 

0.0073 

0.0075 

0.0076 

0.0078 


Additional file 2 — Global sensitivity analysis results for 
lacuna count 

Table S2. Global sensitivity analysis results for lacuna count. 


P 

11 

12 

13 

Variance data 

15.8651 

15.8651 

15.8651 

Variance PGE 

15.5295 

15.9661 

16.7974 

S(Ae) 

0.0075 

0.0074 

0.0070 

S{D) 

0.7120 

0.7130 

0.7187 

S(Aa) 

0.0129 

0.0125 

0.0119 

S((/cell,cell) 

0.0418 

0.0407 

0.0387 

S(Ac,//) 

0.0579 

0.0570 

0.0552 

S(Ac,Aa) 

0.0043 

0.0043 

0.0041 

S( Ac, (/cell, cell) 

0.0144 

0.0145 

0.0143 

SiDM) 

0.0350 

0.0347 

0.0348 

S(//, (/cell,cell) 

0.0533 

0.0521 

0.0501 

S( Ayi,(/cell,cell) 

0.0050 

0.0048 

0.0046 

SiXc.DM) 

0.0289 

0.0315 

0.0331 

S( Ac,//, (/cell, cell) 

0.0222 

0.0232 

0.0236 

S (Ac, A , (/cell,cell) 

0.0118 

0.0131 

0.0142 

S(//,A^, (/cell, cell) 

0.0086 

0.0094 

0.0100 
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Additional file 3 — Global sensitivity analysis results for 
lacuna count capped at maximum of 5 lacunae 

Table S3 Global sensitivity analysis results for lacuna count capped at maxi- 
mum of 5 lacunae._ 


p 

12 

13 

14 

Variance data 

3.4870 

3.4870 

3.4870 

Variance PCE 

3.4346 

3.4821 

3.5668 

S(Ae) 

0.0081 

0.0081 

0.0079 

S{D) 

0.6893 

0.6861 

0.6847 

S(A^) 

0.0078 

0.0077 

0.0077 

S(</cell,cell) 

0.0578 

0.0571 

0.0557 

SiXc.D) 

0.0561 

0.0557 

0.0547 

SiXc.XA) 

0.0084 

0.0085 

0.0086 

S( Ac 5 (/cell, cell) 

0.0394 

0.0409 

0.0429 

S{DM) 

0.0048 

0.0048 

0.0048 

S(//, (/cell,cell) 

0.0489 

0.0485 

0.0476 

S( Ayi,(/cell,cell) 

0.0035 

0.0035 

0.0035 

S(Ac,i/,AA) 

0.0147 

0.0160 

0.0169 

S( Ac,//5 (/cell, cell) 

0.0458 

0.0476 

0.0490 

S(Ac,AA 5 (/cell,cell) 

0.0239 

0.0258 

0.0274 

S(//,A^, (/cell,cell) 

0.0121 

0.0132 

0.0140 
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Additional file 4 — Two-dimensional intensity plots of la¬ 
cuna count with maximum intensity of 5 lacunae 
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Figure SI. Two-dimensional intensity plots of lacuna count with 
maximum intensity of 5 The intensity of the output measure lacuna count, 
mapped to an interval of 0 to 5 as indicated by the color bars, is plotted for 
each two-parameter combination of the parameters the cell rigidity (A, 4 ), cell¬ 
cell adhesion (Jceii.ceii), the diffusion coefficient of the chemoattractant (D), and 
sensitivity of cells to the chemoattractant at cell-matrix interfaces (Ac). 
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